An integer formula for Fibonacci numbers

This code, somewhat surprisingly, generates Fibonacci numbers.

def fib(n):
    return (4 << n*(3+n)) // ((4 << 2*n) - (2 << n) - 1) & ((2 << n) - 1)

In this blog post, I’ll explain where it comes from and how it works.

Before getting to explaining, I’ll give a whirlwind background overview of Fibonacci numbers and how to compute them. If you’re already a maths whiz, you can skip most of the introduction, quickly skim the section “Generating functions” and then read “An integer formula”.

Overview

The Fibonacci numbers are a well-known sequence of numbers: $$1,1,2,3,5,8,13,21,34,55,\dots$$

The $n$th number in the sequence is defined to be the sum of the previous two, or formally by this recurrence relation: $$\begin{array}{rcl}\mathrm{Fib}(0)& =& 1\\ \mathrm{Fib}(1)& =& 1\\ \mathrm{Fib}(n)& =& \mathrm{Fib}(n-1)+\mathrm{Fib}(n-2)\end{array}$$

I’ve chosen to start the sequence at index 0 rather than the more usual 1.

Computing Fibonacci numbers

There’s a few different reasonably well-known ways of computing the sequence. The obvious recursive implementation is slow:

def fib_recursive(n):
    if n < 2: return 1
    return fib_recursive(n - 1) + fib_recursive(n - 2)

An iterative implementation works in $O(n)$ operations:

def fib_iter(n):
    a, b = 1, 1
    for _ in xrange(n):
        a, b = a + b, a
    return b

And a slightly less well-known matrix power implementation works in $O(\log n)$ operations.

def fib_matpow(n):
    m = numpy.matrix('1 1 ; 1 0') ** n
    return m.item(0)

The last method works by considering the a and b in fib_iter as sequences, and noting that $$\left(\begin{array}{c}{a}_{n+1}\\ b_{n+1}\end{array}\right)=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}{a}_n\\ b_n\end{array}\right)$$

From which follows $$\left(\begin{array}{c}{a}_n\\ b_n\end{array}\right)={\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)}^n\left(\begin{array}{c}1\\ 1\end{array}\right)$$

and so if $m={\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)}^n$ then $b_n=m_{11}$ (noting that unlike Python, matrix indexes are usually 1-based).

It’s $O(\log n)$ based on the assumption that numpy’s matrix power does something like exponentation by squaring.

Another method is to find a closed form for the solution of the recurrence relation. This leads to the real-valued formula: $\mathrm{Fib}(n)=({\varphi }^{n+1}-{\psi }^{n+1})/\sqrt{5}$ where $\varphi =(1+\sqrt{5})/2$ and $\psi =(1-\sqrt{5})/2$. The practical flaw in this method is that it requires arbitrary precision real-valued arithmetic, but it works for small $n$.

def fib_phi(n):
    phi = (1 + math.sqrt(5)) / 2.0
    psi = (1 - math.sqrt(5)) / 2.0
    return int((phi ** (n+1) - psi ** (n+1)) / math.sqrt(5))

Generating Functions

A generating function for an arbitrary sequence $a_n$ is the infinite sum ${\mathrm{\Sigma }}_n{a}_n{x}^n$. In the specific case of the Fibonacci numbers, that means ${\mathrm{\Sigma }}_n\mathrm{Fib}(n)x^n$. In words, it’s an infinite power series, with the coefficient of $x^n$ being the $n$th Fibonacci number.

Now, $\mathrm{Fib}(n+2)=\mathrm{Fib}(n+1)+\mathrm{Fib}(n)$

Multiplying by $x^{n+2}$ and summing over all $n$, we get: $${\mathrm{\Sigma }}_n\mathrm{Fib}(n+2)x^{n+2}={\mathrm{\Sigma }}_n\mathrm{Fib}(n+1)x^{n+2}+{\mathrm{\Sigma }}_n\mathrm{Fib}(n)x^{n+2}$$

If we let $F(x)$ to be the generating function of $\mathrm{Fib}$, which is defined to be ${\mathrm{\Sigma }}_n\mathrm{Fib}(n)x^n$ then this equation can be simplified: $$F(x)-x-1=x(F(x)-1)+x^{2}F(x)$$

and simplifying, $$F(x)=xF(x)+x^{2}F(x)+1$$

We can solve this equation for $F$ to get $$F(x)=\frac{1}{1-x-x^2}$$

It’s surprising that we’ve managed to find a small and simple formula which captures all of the Fibonacci numbers, but it’s not yet obvious how we can use it. We’ll get to that in the next section.

A technical aside is that we’re going to want to evaluate $F$ at some values of $x$, and we’d like the power series to converge. We know the Fibonacci numbers grow like ${\varphi }^n$ and that geometric series ${\mathrm{\Sigma }}_n{a}^n$ converge if $|a|<1$, so we know that if $|x|<1/\varphi \simeq 0.618$ then the power series converges.

An integer formula

Now we’re ready to start understanding the Python code.

To get the intuition behind the formula, we’ll evaluate the generating function $F$ at ${10}^{-3}$. $$\begin{array}{rl}& F(x)=\frac{1}{1-x-x^2}\\ & F({10}^{-3})=\frac{1}{1-{10}^{-3}-{10}^{-6}}\\ & =1.001002003005008013021034055089144233377610988599588187\dots \end{array}$$

Interestingly, we can see some Fibonacci numbers in this decimal expansion: $1,1,2,3,5,8,13,21,34,55,89$. That seems magical and surprising, but it’s because $F({10}^{-3})=\mathrm{Fib}(0)+\mathrm{Fib}(1)/{10}^3+\mathrm{Fib}(2)/{10}^6+\mathrm{Fib}(3)/{10}^9+\dots$.

In this example, the Fibonacci numbers are spaced out at multiples of $1/1000$, which means once they start getting bigger that 1000 they’ll start interfering with their neighbours. We can see that starting at 988 in the computation of $F({10}^{-3})$ above: the correct Fibonacci number is 987, but there’s a 1 overflowed from the next number in the sequence causing an off-by-one error. This breaks the pattern from then on.

But, for any value of $n$, we can arrange for the negative power of 10 to be large enough that overflows don’t disturb the $n$th Fibonacci. For now, we’ll just assume that there’s some $k$ which makes ${10}^{-k}$ sufficient, and we’ll come back to picking it later.

Also, since we’d like to use integer maths (because it’s easier to code), let’s multiply by ${10}^{kn}$, which also puts the $n$th Fibonacci number just to the left of the decimal point, and simplify the equation. $${10}^{kn}F({10}^{-k})=\frac{{10}^{kn}}{1-{10}^{-k}-{10}^{-2k}}=\frac{{10}^{kn+2k}}{{10}^{2k}-{10}^k-1}$$

If we take this result modulo ${10}^k$, we’ll get the $n$th Fibonacci number (again, assuming we’ve picked $k$ large enough).

Before proceeding, let’s switch to base 2 rather than base 10, which changes nothing but will make it easier to program. $$2^{kn}F(2^{-k})=\frac{2^{k(n+2)}}{2^{2k}-2^k-1}$$

Now all that’s left is to pick a value of $k$ large enough so that $\mathrm{Fib}(n+1)<2^k$. We know that the Fibonacci numbers grow like ${\varphi }^n$, and $\varphi <2$, so $k=n+1$ is safe.

So! Putting that together: $$\begin{array}{rl}\mathrm{Fib}(n)& \equiv 2^{(n+1)n}F(2^{-(n+1)})\mathrm{mod}{2}^{n+1}\\ & \equiv \frac{2^{(n+1)(n+2)}}{(2^{n+1}{)}^2-2^{n+1}-1}\mathrm{mod}{2}^{n+1}\\ & \equiv \frac{2^{(n+1)(n+2)}}{2^{2n+2}-2^{n+1}-1}\mathrm{mod}{2}^{n+1}\end{array}$$

If we use left-shift notation that’s available in python, where $amathrm«k=a\cdot 2^k$ then we can write this as: $$\mathrm{Fib}(n)\equiv \frac{4<<n(3+n)}{(4<<2n)-(2<<n)-1}\mathrm{mod}(2<<n)$$

Observing that $\mathrm{mod}(2«n)$ can be expressed as the bitwise and (&) of $(2«n)-1$, we reconstruct our original Python program:

def fib(n):
    return (4 << n*(3+n)) // ((4 << 2*n) - (2 << n) - 1) & ((2 << n) - 1)

Although it’s curious to find a non-iterative, closed-form solution, this isn’t a practical method at all. We’re doing integer arithmetic with integers of size $O(n^2)$ bits, and in fact, before performing the final bit-wise and, we’ve got an integer that is the first $n$ Fibonacci numbers concatenated together!